#include <vector>
#include <boost/numeric/interval.hpp>

using namespace boost::numeric;

struct point
{
   double x;
   double y;
};

std::pair<double, double> sum(double a, double b)
{
   std::pair<double, double> res;
   double x, y;
   x = a + b;
   double b_v = x - a;
   double a_v = x - b_v;
   double b_r = b - b_v;
   double a_r = a - a_v;
   y = a_r + b_r;
   res.first = x;
   res.second = y;
   return res;
}

std::pair<double, double> split(double a, int s)
{
   std::pair<double, double> res;
   double c = ((1LL << s) + 1LL) * a;
   double a_big = c - a;
   res.first = c - a_big;
   res.second = a - (c - a_big);
   return res;
}

std::pair<double, double> product(double a, double b)
{
   std::pair<double, double> res;
   int p = std::numeric_limits<double>::digits;
   double x, y;
   x = a * b;
   int p2 = p / 2;
   if(p2 * 2 != p)
      p2++;

   double a_hi, a_lo, b_hi, b_lo;

   a_hi = split(a, p2).first;
   a_lo = split(a, p2).second;
   b_hi = split(b, p2).first;
   b_lo = split(b, p2).second;

   double err1, err2, err3;
   err1 = x - a_hi * b_hi;
   err2 = err1 - a_lo * b_hi;
   err3 = err2 - a_hi * b_lo;
   y = a_lo * b_lo - err3;
   res.first = x;
   res.second = y;

   return res;
}


void grow_expansion(std::vector<double> &e, double b)
{
   double q = b;
   for(int i = 0; i < e.size(); i++)
   {
      std::pair<double, double> s = sum(q, e[i]);
      e[i] = s.second;
      q = s.first;
   }
   e.push_back(q);
}


std::vector<double> expansion_sum(const std::vector<double> &e, const std::vector<double> &f)
{
   std::vector<double> res = e;
   for(int i = 0; i < f.size(); i++)
   {
      double q = f[i];
      for(int j = i; j < res.size(); j++)
      {
         std::pair<double, double> s = sum(q, res[j]);
         res[j] = s.second;
         q = s.first;
      }
      res.push_back(q);
   }

   return res;
}

int left_turn(point a, point b, point c)
{
   double eps_m = std::numeric_limits<double>::epsilon();
   double v = std::abs((b.x - a.x) * (c.y - a.y)) + std::abs((b.y - a.y) * (c.x - a.x));
   double s = (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
   double eps = 4 * eps_m * v;
   if(s > eps)
      return 1;
   if(s < -eps)
      return -1;

        //typedef interval<double> id;

        //id res = (id(b.x) - id(a.x)) * (id(c.y) - id(a.y)) -
        //  (id(b.y) - id(a.y)) * (id(c.x) - id(a.x));

        //if(res.lower() > 0)
        //  return 1;
        //if(res.upper() < 0)
        //  return -1;

   std::vector<double> e1, e2, e3;
   std::pair<double, double> p;

   std::pair<double, double> v_arr[6];

   v_arr[0] = product(b.x, c.y);
   v_arr[1] = product(b.x, -a.y);
   v_arr[2] = product(-a.x, c.y);

   v_arr[3] = product(-b.y, c.x);
   v_arr[4] = product(b.y, a.x);
   v_arr[5] = product(a.y, c.x);

   std::vector<double> exp;
   for(int i = 0; i < 6; i++)
   {
      grow_expansion(exp, v_arr[i].first);
      grow_expansion(exp, v_arr[i].second);
   }

   for(int i = 11; i >= 0; i--)
   {
      if(exp[i] > 0)
         return 1;
      if(exp[i] < 0)
         return -1;
   }
   return 0;
}


bool bounding_box(double a, double b, double c, double d)
{
   if(a > b)
   {
      double tmp = a;
      a = b;
      b = tmp;
   }

   if(c > d)
   {
      double tmp = c;
      c = d;
      d = tmp;
   }

   if(a > d || c > b)
      return false;
   return true;
}

bool intersection(point a, point b, point c, point d)
{
   return bounding_box(a.x, b.x, c.x, d.x) && bounding_box(a.y, b.y, c.y, d.y) &&
          left_turn(a, b, c) * left_turn(a, b, d) <= 0 && left_turn(c, d, a) * left_turn(c, d, b) <= 0;
}

int horizontalIntersection(point a, point b, double line)
{
   if(bounding_box(a.y, b.y, line, line))
   {
      if(a.y == line || b.y == line)
      {
         return 0;
      }
      return 1;
   }
   return -1;
}